Nonlinear resonant tunneling of Bose-Einstein condensates in tilted optical lattices 
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We study the tunneling decay of a Bose-Einstein condensate out of tilted optical lattices within 
the mean-field approximation. We introduce a novel method to calculate also excited resonance 
eigenstates of the Gross-Pitaevskii equation, based on a grid relaxation procedure with complex ab- 
sorbing potentials. This algorithm works efficiently in a wide range of parameters where established 
methods fail. It allows us to study the effects of the nonlinearity in detail in the regime of resonant 
tunneling, where the decay rate is enhanced by resonant coupling to excited unstable states. 
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I. INTRODUCTION 

The dynamics of a quantum particle in a periodic po- 
tential subject to an external force is one of the central 
problems in solid state physics. In the field free case all 
eigenstates are delocalized over the lattices, leading to 
transport [l|, Qj . The application of a constant force leads 
to a localization of the eigenstates such that transport is 
suppressed contrary to our intuition [J-Q- Instead, the 
quantum particle performs the celebrated Bloch oscilla- 
tions, and eventually decays by repeated Zener tunneling 
to higher Bloch bands |7H17|. The most detailed stud- 
ies of Bloch oscillations and decay have been carried out 
with ultracold atoms trapped in optical lattices. These 
systems are particularly appealing, because the dynamics 
of the atoms can be recorded in situ and all parameters 
can be tuned precisely over a wide range. The external 
force can be induced by gravity [7|, magnet ic gradient 
fields [3] or by accelerating the lattice |9l-fl3{. Decay in 
strong fields manifests itself in the pulsed output of coher- 
ent matter waves. The dynamics is even more interesting 
when the atoms undergo Bose-Einstein condensation and 
interactions have to be taken into account. For low tem- 
perature and high densities, the dynamics of the atoms 
can be described by the celebrated Gross-Pitaevskii equa- 
tion (GPE) with astonishing accuracy fls] . In this treat- 
ment, interactions are incorporated by a nonlinear mean- 
field potential, which is proportional to the condensate 
density. The nonlinearity of the equation alters the dy- 
namics and in particular the decay substantially. In- 
teractions can lead to a damping of Bloch oscillations 
|19l] . asymmetric Landau- Zener tunneling [ifl, [20) [111 , or 
a bistability of resonance curves |22l424| . 

Here we study the resonance eigenstates of the GPE 
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with a periodic potential V{x -\- d) = V{x) and a static 



force F > 0, which is known as a Wannier- Stark (WS) 
potential. The imaginary part F of the eigenenergy gives 
the decay rate of the condensate. In the following we 
use scaled units with h — m — \ and we consider a 
cosine potential V{x) — cos(a;) unless otherwise stated. 
A comprehensive review of the localized eigenstates, the 
WS resonances, can be found in [16i] . 



In this article we introduce a new algorithm for the 
computation of nonlinear resonance states based on a 
grid relaxation method with a complex absorbing poten- 
tial (CAP). This algorithm converges in a wide parameter 
range and is applicable even to situations of many degen- 
erate energy levels, such as the WS system at resonance 
condition (see below). It is thus capable to describe gen- 
uine nonlinear phenomena such as bistability, which pose 
a major difficulty to other methods as for instance non- 
linear complex scaling (CS) |25l - [29| . In addition, it is 
more efficient and easier to implement and, unlike pre- 
vious methods, is not restricted to ground state calcula- 
tions but can also compute excited states. Note that our 
approach differs from the CAP method used in [2g, |30| 
because the latter does not use a grid relaxation but relies 
on a basis set expansion. Though such expansions work 
well for simple single well potentials, they cannot easily 
handle complicated problems like the Wannier-Stark sys- 
tem studied in the present paper, which requires the use 
of as much as 500 basis states even in the linear (nonin- 
teracting) case [31|. Our method is applied to study the 
decay of a Bose-Einstein condensate in the strongly non- 
linear regime. Nonlinear effects are crucial in the regime 
of resonantly enhanced tunneling (RET). In this case a 
metastable WS resonance becomes energetically degener- 
ate with an excited, less stable state, which can increase 
the decay rate by orders of magnitude. This phenomenon 
is most pronounced in deep optical lattices and has been 
studied systematically for the linear case in ^15, il6| . The 
nonlinearity shifts the resonance and eventually bends 
the resonance peak leading to a bistable behavior. 



II. COMPUTATIONAL METHOD 

Linear WS resonances can be efficiently calculated with 
the truncated shift operator technique introduced in [31| . 
In the nonlinear case, the method of CS has been applied 
[25l [27h29| . Though satisfactory from a conceptual point 
of view, this method has several drawbacks. The im- 
plementation is complicated as it requires switching be- 
tween different basis sets as well as different time propa- 
gation methods. Furthermore, the calculation of excited 
states is highly non-trivial, as the method relies on an 
imaginary time propagation, and the convergence is quite 
slow, especially for weak fields and close to energetic de- 
generacies as present in the RET condition ^, J8] . As 
an alternative, we propose a method based on complex 
absorbing potentials (CAP) performed on a finite grid 
[a;_,a;+] in real space. We assume that the resonance 
wave function is mainly localized in the interval [a:^,^;^] 
with X- < Xi < Xr < a;+ and fix the normalization as 
/ '^ \ip{x)\'^dx = 1. For X — > — cxD, we apply a CAP of the 
type 



VcAP oc 







X < X- 
X > X- 



(2) 



which only modifies the wave function in the vicinity of 
the grid boundary x- making it square integrable. For 
X —^ -|-oo, the wavefunction rapidly converges to zero, so 
that no further CAP is needed on this side. The bound- 
ary conditions for the wave function read 



V'(x_)=0, ^(x+)=0, ^'{x+)=C, 



(3) 



where the last condition is used to control the normaliza- 
tion. The algorithm starts from the linear case g = 0, for 
which all WS resonances can be computed efficiently [3l| . 
Nonlinear WS resonances in different bands are calcu- 
lated by choosing a different initial guess. The nonlinear- 
ity is then increased gradually, using the previous result 
as initial guess for a standard boundary value problem 
(BVP) solver, e.g. the Matlab -function bvp4c. Apply- 
ing the BVP solver changes the normalization of ip, such 
that the parameter C has to be adjusted according to 



C-^C 



1-0(3;) I dx 
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(4) 



This is repeated until the normalization converges to 
unity. Then the nonlinearity is increased by one step. To 
demonstrate the validity of this algorithm we compare 
the calculated decay rates for a cosine potential for sev- 
eral parameters to complex scaling results, which them- 
selves were tested against a direct time propagation in 
Ref. [231. The values summarized in Tab. U show an ex- 
cellent agreement over the entire parameter range. Resid- 
ual numerical errors are very small; they can mainly be 
attributed to the limited computation time for the CS 
method and refiections of the matter wave at the CAP. 
For a further discussion of CAPs in the simulation of few 
boson systems, see [32] and references therein. 



9 


F 


Pes 


PcAP 







0.5 


1.941 X 10-2 


1.941 


xlO-2 


0.1 


0.5 


2.180 X 10-2 


2.180 


X 10-2 





0.25 


7.2 X 10"* 


7.104 


xlO-* 


0.1 


0.25 


8.4 X 10-* 


8.346 


xlO-* 


0.2 


0.25 


9.7 X 10-* 


9.688 


xlO-* 


0.25 


0.25 


1.04 X 10-^ 


1.041 


xlO-3 


0.5 


0.25 


1.48 X 10-^ 


1.476 


X 10-^ 


0.2 


0.15 


2.9 X 10-^ 


2.832 


X 10-^ 


0.2 


0.13125 


5.7 X 10-^ 


5.600 


X lO-'' 



TABLE I: Decay rates T for the most stable resonance of the 
potential V{x) = cos(a;), taken from Ref. [23| (CS method) 
and computed by the CAP grid relaxation method. Particu- 
larly for small decay rates the new CAP method proves more 
efficient than the CF technique. 



III. RESONANTLY ENHANCED TUNNELING 

We use the CAP method to investigate how a nonlinear 
interaction affects the decay of a BEC in a tilted optical 
lattice. In the weakly interacting regime, the scaling of 
the decay rate with the field strength is given by the cel- 
ebrated Landau-Zener formula T{F) « i^exp(7rAi?^/i^), 
where AE is the energy gap between the Bloch bands of 
the periodic potential 0, [iJ] and the field strength F de- 
termines the oscillation frequency in the bands [IJ, [la] . 
Major differences arise in the regime of RET. In this case 
an eigenstate localized mainly in one of the wells of the 
potential becomes energetically degenerate with an ex- 
ited state in another well, which can increase the decay 
rate by orders of magnitude flfjl]. In the following, we fo- 
cus on the experimentally studied regime [Uj, [l2| , where 
already a modest nonlinearity strongly affects the decay 
of the condensate [ill. [T^. [TtJ . 

RET is illustrated in Fig. [1] (a,b) for the linear case 
(7 = 0, showing the decay rate P and the chemical po- 
tential /i of the two most stable resonances as a function 
of F. RET is observed at 1/F w 7.5, where the two 
energy levels /x(i^) cross. The resonant coupling to the 
excited states leads to a pronounced RET peak of the 
decay rate for the most stable resonance. Coincidentally, 
a pronounced dip is observed for the first excited reso- 
nance, which is stabilized by the coupling to the most 
stable resonance [16|. The influence of a small nonlin- 
earity is illustrated in Fig. [1] (c). Three main effects are 
observed: a shift of the resonance peaks, an increase (de- 
crease) of the peak decay rate in the ground state for 
g > {g < 0) and a deformation of the peak shape. 

The shift and the deformation can be qualitatively un- 
derstood by a perturbative approach ^27]. To first order, 
this predicts a shift of the real part of the eigenenergy. 



AfM{g) « g / \i,g\'\i,g^ordx « g / l^g=o|"dx, (5) 

'^Xi J Xl 

which corresponds to a shift of the field strength accord- 





FIG. 1: (Color online) Resonantly enhanced tunneling (RET) 
of (non)-linear WS resonances, (a) Energies and (b) de- 
cay rates of the two most stable WS resonances in a cosine- 
potential as a function of the inverse field strength 1/F. (c) 
Shift of the RET peaks due to the nonlinear interaction of 
a BEC for g = -h0.02 (o), g = -0.02 (o) and p = (- -). 
Numerical results (symbols) are compared to a perturbative 
calculation (solid lines) according to Eqs. ((SJ and ((TJ. 



ing to 



^F{g) « ±^^l{g)/{2^T). 



(6) 



Here, the minus sign holds for the ground and the plus 
sign for the excited band. The nonlinear decay rate is 
then approximately given by 



Tg{F)^T^{F- 



^F{9))- 



(7) 



The shift is further investigated in Fig. [5] (b), where 
the decay rate as well as the peak position is plotted vs. 
the interaction strength over a wide parameter range. 
The perturbative calculation ([S|) predicts that the peak 
position Ft-cs is shifted with a slope dF,-cs/dg = 0.059 
for small values of g, which is plotted as a green line 
in Fig. [5] (b) . This deviates from the numerically exact 
results already for small values of g, for which a linear fit 
yields a smaller slope of dF^e^/dg = 0.051. In agreement 
with [27| we thus find that first order perturbation theory 
is not sufficient to describe the shift of the RET peaks 
quantitatively. Noticeably, the RET peak and the dip of 
the decay rate for the first excited resonance always shift 
into opposite directions, as shown in Fig.[l](c). 

The change in the maximum decay rate is not predicted 
by perturbation theory, but easily explained phenomeno- 
logically. It is a direct consequence of the interaction 
as repulsion between the particles in general leads to a 
destabilization whereas attraction leads to a stabilization 
of both resonances and bound states (see [l3j and refer- 
ences therein). This is further illustrated in Fig. [5] (c), 



FIG. 2: (Color online) (a) Colormap plot of the decay rates of 
the most stable WS resonances in a cosine-potential vs. the 
field strength F and the interaction strength g in the vicinity 
of the first order RET peak, (b) Position and (c) height of 
the RET peak vs. the interaction strength g. 



where the peak decay rate of the most stable resonance 
is plotted as a function oi g over a wide parameter range. 
Similar effects have been investigated for several other 
model potentials [il, [H, [H . 

Figure[T](c) furthermore shows that the RET peaks be- 
come asymmetric for g ^ Q. For a repulsive (attractive) 
nonlinearity, the peak bends to higher (lower) values of 
F . If the nonlinearity is increased above a critical value 
5cr, the peaks bend over and a bistable behavior emerges 
as shown in Fig. [3l For the given parameters of the op- 
tical lattice, the critical nonlinearity has been found at 
0.5 < 5cr < 0.6. The detailed shape of a bistable RET 
peak is plotted in the upper panel of Fig. [31 which also 
indicates how WS states are calculated numerically in the 
bistable regime: We have started with a small value of 
F which was then gradually increased, using every result 
as initial guess for the next calculation (cf. also |27ll28|). 
After reaching a final, large value of the field strength, 
the procedure was reversed and F was decreased back to 
the initial value. Within the regime of bistability forward 
and backward sweep yield the upper and lower branches 
of the peak, respectively. The intermediate branch is gen- 
erally difficult to compute as it is dynamically unstable. 

The bending of the RET peak and the emergence of a 
bistability can be understood qualitatively by the pertur- 
bative approach introduced above. A common WS state 
in a deep optical lattice is strongly localized in a single 
potential well so that its chemical potential is strongly 
changed according to Eq. d?]). In comparison, the state 
corresponding to the maximum of the RET peak is de- 
localized because of the energetical degeneracy with an 
excited state in another well. Therefore its chemical po- 
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FIG. 3: (Color online) Upper panel: Bistability of the RET 
peak for strong repulsive interactions {g = 0.8). The decay 
rate was calculated for a forward sweep (blue asterisk) and a 
backward sweep (red circles). A spline interpolation (dashed 
line) is included to guide the eye. The solid line shows the 
linear {g = 0) peak shape for comparison. Lower panel: Emer- 
gence of the bistable behavior: Decay rate as a function of the 
field strength F for difi'erent values of the nonlinearity g. 



tential is afTected rather weakly and according to Eq. ^ 
also the change of the peak position AF is small (cf. also 
[111 [l3|). With increasing nonlinearity, the edges of a 
RET peak shift to smaller values of the field strength, 
while the maximum falls behind. The whole peak bends 
to the right and finally becomes bistable. 

Bifurcations of nonlinear stationary states have been 
previously observed in several different contexts: The 
simplest example is the nonlinear two-mode system. As 
observed for the WS state, new stationary states emerge 
first in the vicinity of the resonance, i.e. the when the 
chemical potential of the two coupled states coincides 
(see [20, l22j and references therein) . The emergence of 
bistability has also been analyzed for the transmission 
coefficient in the context of nonlinear RET through one- 
dimensional potential barriers 22-24]. However, in this 
case states corresponding to the transmission maximum 
are localized strongest. Thus the resonance peaks bend 
into the same direction as they are shifted, which is op- 
posite from the behavior of the WS RET peaks shown in 
Fig.O 
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FIG. 4: (Color online) Chemical potential /i and decay rates F 
for (non) linear WS-resonances in a bichromatic optical lattice 
for (5 = 1 and (a) (jj = 3.4, (b) <j) = 3.2, (c) (jj = 3.292 and 
g = (dashed black line), g — -1-0.02 (thick blue line) and 
g = —0.02 (dash-dotted red line). 



IV. BEYOND THE RET-REGIME 

A new regime of RET can be explored in bichromatic 
optical lattices, 

1 A 

V{x) = -cos{x) + -cos{2x + (t)). (8) 

These potentials can be realized experiment ally by su- 
perimposing two incoherent optical lattices [3J, [S^, or 
by combining optical potentials based on virtual two- 
photon and four-photon processes [3a,[3g. It can be used 
to study Landau-Zener tunneling between different bands 
and the interplay of tunneling and Bloch oscillations [331 ■ 
Depending on the relative phase and the relative 
strength of the two lattices, WS resonances in tilted 
bichromatic optical lattices show two remarkably differ- 
ent types of level crossing scenarios [33 - either the real 
parts fi (type-I) or the imaginary parts F (type-II) of the 
eigenenergies cross. A full degeneracy of both ^ and F 
occurs only for isolated points in parameter space. Ex- 
amples are shown in Fig.[3|for ^ = 1 and different relative 
phases 4> of the two lattices. A familiar type-II crossing is 
observed for (p = 3.4. The real parts of the eigenenergies 
(/i) cross, while the imaginary parts F anti-cross, leading 
to the familiar RET-peaks of the decay rates. Changing 
the phase slightly to (j) — 3.2, one finds a type-I cross- 
ing. The decay rates of the two most stable resonances 
now cross, while the real parts anti-cross. An exceptional 
point [38j, where both real and imaginary part are fully 
degenerate, is found for (j> = 3.292, as shown in Fig. [3| 



(c). However, the degeneracy is lifted as soon as the 
atoms start to interact. A weak repulsive nonlinearity 
g — +0.02, turns the exceptional crossing into an or- 
dinary type-I crossing, while an attractive nonlinearity 
g = —0.02 favors a type-II crossing. This change of peak 
shape can have dramatic effects on the dynamics of a 
Bose-Einstein condensate, in particular when experimen- 
tal parameters are adiabatically varied (see, e.g. [38]). 



CONCLUSIONS 



nant tunneling, where the decay rate can be enhanced 
by orders of magnitude by resonant coupling to unsta- 
ble excited states. The interactions shift and bend the 
resonance peaks and eventually lead to a bistable peak 
shape. Even more interesting effects can be studied in 
tilted bichromatic lattices, where different types of level 
crossing scenarios emerge when the lattice parameters are 
tuned. These effects will be studied in detail in a future 
publication. 



Bose-Einstein condensates in tilted optical lattices are 
ideal to study the decay of interacting open quantum sys- 
tems. Experimentally the parameters can be tuned over 
a wide range and the dynamics can be recorded in situ. 
Here we presented an efficient method to calculate the 
decay rate in the mean-field regime also in the presence 
of degeneracies which also, unlike previous methods, is 
not restricted to ground state calculations. The effects 
of the nonlinearity are strongest in the regime of reso- 
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